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Abstract 

The numerical instability observed in the Electromagnetic-Particle-in-cell 
(EM-PIC) simulations with a plasma drifting with relativistic velocities is 
studied using both theory and computer simulations. We derive the numer- 

C\| ical dispersion relation for a cold plasma drifting with a relativistic velocity 

and find an instability attributed to the coupling between the beam modes of 
the drifting plasma and the electromagnetic modes in the system. The char- 
acteristic pattern of the instability in Fourier space for various simulation se- 
tups and Maxwell Equation solvers are explored by solving the corresponding 
numerical dispersion relations. Furthermore, based upon these characteristic 

CN patterns we derive an asymptotic expression for the instability growth rate. 

The asymptotic expression greatly speeds up the calculation of instability 
growth rate and makes the parameter scan for minimal growth rate feasible 

^ even for full three dimensions. The results are compared against simulation 

results and good agreement is found. These results can be used as a guide 
to develop possible approaches to mitigate the instability. We examine the 
use of a spectral solver and show that such a solver when combined with a 
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low pass filter with a cutoff value of \k\ essentially eliminates the instability 
while not modifying modes of physical interest. The use of spectral solver 
also provides minimal errors to electromagnetic modes in the lowest Brillouin 
zones. 

Keywords: Particle-in-cell, plasma simulation, relativistic plasma drift, 
numerical dispersion relation, numerical instability, numerical Cherenkov 
radiation 



1. Introduction 

The effect of finite grid size and finite time step in electromagnetic sim- 
ulations using the particle-in-cell method has been extensively studied both 
theoretically and in simulations [H |2j |3j 0]. Understanding these effects 
are crucial when one tries to separate numerical artifacts from real plasma 
phenomena. In recent years simulations with plasmas drifting with relativis- 
tic speeds have been conducted extensively, owing to using Lorentz boosted 
frames for the study of Laser Wakefield Acceleration (LWFA) [7J IHl El [TO] 
and for the study of relativistic collisionless shocks. These simulations have 
revealed a numerical instability which only occurs in multi-dimensions and 
which limits the range of parameters that can be explored in Lorentz boosted 
frame simulations of LWFA and relativistic shocks (TTj [T2]. Godfrey stud- 
ied the numerical instability induced by the a plasma drifting in ID [3] and 
in mult i- dimensions [4j. These analysis did not include relativistic mass ef- 
fects in the plasma response (and therefore did not include relativistic drifts) 
and were applied to a code that solved the scalar and vector potential using 
the Coulomb gauge. However, most present day codes solve directly for the 
electric and magnetic fields and use a rigorous charge conserving current de- 
position. Recent results indicate that there is only an instability in 2D and 

3D judo]. 

To better understand and with the hope of mitigating the observed insta- 
bility, we present here an analysis of the numerical instability for a plasma 
drifting relativistically in multi-dimensions in which the electric field E and 
magnetic field B are solved for directly and where only a current deposit is 
included. We follow the basic method and notation of |lj and concentrate on 
situations where the plasma is cold but is drifting near the speed of light. The 
dispersion relation can be applied for various Maxwell field solvers and it in- 
cludes finite size particle and aliasing effects. The dispersion relation reduces 
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to purely longitudinal plasma waves and purely transverse light waves in a 
drifting plasma in the appropriate limits. The dispersion relation predicts 
the growth rates and range (pattern) of unstable modes in Fourier space. By 
comparing the theoretical predictions against the simulation results using the 
EM-PIC FDTD code OSIRIS |T3] and the UCLA PIC Framework [H] which 
is based on spectral solver, we conclude that the observed instability is in- 
deed induced by the relativistic plasma drift (and not due to under-resolved 
backscattered radiation). It is found that the unstable modes lie near the 
intersection between beam modes and transverse EM modes. We use this 
fact to develop asymptotic expressions for the instability growth rate. These 
observations can then be used as a guide for selecting alternative Maxwell 
Equation solvers and smoothing schemes to mitigate the instability. Specif- 
ically, we show that a spectral solver together with a cutoff filter in k space 
can eliminate the instability. 

The remainder of this paper is organized as follows. We first derive the 
multi-dimensional numerical dispersion relation in section [2] In section [3j we 
use the 2D dispersion relation obtained in section [2] to study the numerical 
instability induced by relativistic plasma drift. The theoretical results are 
compared against the simulation results, and good agreement is found. After 
including a minor correction in an earlier version our dispersion relation now 
predicts the time steps for which the minimal instability growth rate was 
observed after reading Ref. [5J. We then use these results as a guide to 
discuss the methods for mitigating this instability. It is found that by using 
a spectral solver to advance the EM field in Fourier space (which supports 
light waves with phase velocity greater than plasma drifting velocity), we 
can obtain desirable patterns of the instability which makes the mitigation 
of it more convenient. To better explore the instability in 3D scenario, we 
developed an asymptotic form of the instability growth rate in section |5j and 
used it to explain the observation that minimal growth rates occur under 
certain time steps, as reported in [8]. Conclusions and direction for future 
work are presented in section [6] The detailed form of the field interpolation 
functions, as well as the finite difference operator used in this paper are listed 
in Appendix A 
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2. Numerical dispersion relation for cold plasma drift 

2.1. Derivation of dispersion relation 

We mainly follow the notation in Ref. [I] to derive the numerical disper- 
sion relation for a cold plasma drifting with relativistic speeds. We note that 
the multi-dimensional analysis in [2] solves for (f> and A and is not valid for 
relativistic drifts, and the ID analysis in |3] predicts growth in ID. On the 
other hand, our analysis includes relativistic mass effects, is valid in multi- 
dimensions, and it predicts instability only in multi-dimensions (in agreement 
with our simulations). Since most EM-PIC codes now in use solve for the 
electric field E and magnetic field B directly (with finite difference or spec- 
tral solvers), we derive a numerical dispersion relation directly using these 
two quantities. Gaussian units will be used; in addition, particle mass and 
velocity will be normalized to electron mass and the speed of light. 

For a multi-dimensional simulation setup in Cartesian coordinates, the 
EM field interpolated on a particle can be expressed as 

E{t, x) = oE(t, m, x, n)E v 



j in. n 



B(t,x) = S B (t,m,x,n)B mi jt 

where m is the time index and n is the grid index; S is the interpolation 
dyadic used to obtain the appropriate field at x and t = mAt; E m ^ and B m ^ 
stands for the electromagnetic forces at time grid index m and space grid 
index n. For momentum conserving field interpolation Se and Sb are equal 
and are scalar functions times the unit dyadic while for energy conserving field 
interpolation and Sb are not eq ual in each dir ection (the interpolation 



dyadics for E, B, and j are given in Appendix A). The momentum change 



of the particle is related to the change in the distribution function of the 
plasma by the linearized Vlasov equation 

d -J(t,x,p) + V - ■ ^f(t,x,p) + q\E(t,x) + V - x B(t,x)\ ■ ^ = 



at 7 ox i 7 j op 

where p is the particle momentum, and 7 is the particle Lorentz factor. After 
Fourier transforming, the Vlasov Equation becomes 



f(u,k,p) = - iq !^E(u,k)E(u,k) + ^x {S~B(cj,k)B(cu,k)}^ ■ ^("-k 



7' 
(1) 
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Note that E and B are defined at the discrete grid position and discrete 
time step, so its Fourier transform in (w, k) is periodic, i.e., 

E(u, k) = E(u', k') B(u, k) = B(u', k') (2) 

where 

tu = to + fxujg ujg = — n = 0, ±1, ±2, . . . 

k'i = k + Vikgi k gi = ^— Vi = 0,±1,±2, ... (3) 

Note that when the EM field are staggered (such as on a Yee mesh), there 
is an addition (— 1)5-* "* term in E(u',k') and B(u>',k'), where % is summed 
over the directions for which the EM fields are half-grid offset [5| ; we absorb 
these additional terms into g aad £ to keep Eq. # eorreet (see^fa 
At. 



Replacing (w, k) with (w', k') in Eq. (jlj), and using Eq. (|3j), we obtain 

/(a/, k',p) = fc) + - x {S(u/, fe)} j • ^(u/ - fe' • 

I 7 J op 7 

(4) 

The current density j due to the movement of the particles can be ex- 
pressed as 

Sj (x' — x)—f({m + 1/2} At, x* ,p)dx'dp 

^ y 

where Sj (x' — x) is the dyadic for the current deposit. After Fourier trans- 
forming we obtain 

^ y 

j(u,k) = qJ2(- 1 )' 1 J SA ~ k ' W f(u'ik',p)dp (5) 

We can now proceed in the normal way to obtain a dispersion relation. 
We start from Faraday's and Ampere's Law, 

Vx£ = -aF 
„ - , -. 

V X 5 = — + 47TJ 

ot 
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which upon Fourier transforming gives, 

[k\ E x E = [lo]B (6) 
[k] B x B = -[u]E-4mj (7) 

[k] E and [k] B are the finite difference operator of the corresponding Maxwell 
solver schemes for E and B fields. We follow the notation in Ref. [3], and 
use [•] exclusively to indicate the finite difference operator. Applying [fc]gX 
to both sides of Eq. (jH), we end up with the coupled wave equation for E 
and j, 

([uj] 2 - [k] E ■ [k] B + [k] E [k] B )E = -Am[uj}j (8) 
Using Eq. <Q and (|5]), we could obtain 

(H 2 -[fc] £ -[fc] B + fi £ fi B )E = -4«g^(-ir[a;] J S ^~ k > /(a/ ,k> ,p)dp 

(9) 



and if we normalize the distribution function such that fo = no/oS use the 
definition of plasma frequency 

u 2 p = 4nq 2 n (10) 

and use the expression for the distribution function in Eq. Q, we finally 
obtain 

([uj] 2 -[k} E -[k} B + [k] E [k] B )E 

.2 



£(_1) J f S A *WP {[ u] F E (u/,Z)E + * x {g(u/, x E)}\ 

(11) 

which is a generalized dispersion relation for a plasma of finite size particles 
drifting on a grid. We note that the use of additional smoothers and filters 
can be incorporated into the dispersion relation by adding additional SsM(k') 
terms outside the summation over Brillouin zones (essentially it multiplies 
the uj 2 term). 



6 



2.2. Elements of dispersion relation tensor 

We next examine the dispersion relation in the limit of a cold plasma 
including the possibility that the drift is near the speed of light. 

2.2.1. 3D case 

Note that V for the fields and the current has only three diagonal el- 



ements Si, S 2 , S 3 in each case. In 3D, we can expand Eq. (11) explicitly 

as 

\[u} 2 - [k] E1 [k] m - [k] E2 [k] B2 - [k} E3 [k] B3 )E 1 + [k\ E1 [k\ Bl E 1 + [k} E1 [k} B2 E 2 + [k} E1 [k] B3 E 3 
([a;] 2 - [k} E1 [k} B1 - [k] E2 [k] B2 - [k] E3 [k] B3 )E 2 + \k\B2\k\B\Ex + [k] E2 [k] B2 E 2 + [k] E2 [k] B3 E 3 
^[cu] 2 - [k) E1 [k]m ~ [k] E2 [k] B2 - [k] E3 [k] B3 )E 3 + [/cj^^Bi^i + [k]m[k] B2 E 2 + [k] E3 [k] B3 E 3/ 

d Pl d P2dP3 ( ( ^ J d d f T 

- - IUJ> ~ k[pi ~ k ' 2P2 ~ ^3) \s% 2 ,) [dfl/dZ 

^[u]S E1 Ei +p 2 S B3 ([k] E1 E 2 - \k\B2E\) + PzS B2 {[k] E iE 3 - [k] E3 E x ) 
^}S E2 E 2 +p 3 S B i([k} E2 E 3 - \k} E3 E 2 ) +PiS B3 {[k] E2 E 1 - [k} E1 E 2 ) 
^[u]S E3 E 3 +PiS B2 {[k} E3 E 1 - \k} El E 3 ) + p 2 S B1 {[k] E3 E 2 - \k} E2 E 3 ) 



This can be rewritten as 



4-» 



en ei2 ei3\ (E\ 

£21 £22 C23 I I E 2 

v£3i £32 £33/ \Ez. 



(12) 



(13) 



where we note that e is not the dielectric tensor. In addition, we are most 
interested in a cold plasma that is drifting. For such a case, the unperturbed 
distribution function is given by 



fZ = 5( Pl -po)5(p 2 )5(p 3 ) 



(14) 



where po = jvq, and v is the drifting velocity of the plasma. Substituting the 



above form for f , Eq. (14), into Eq. (12), and carrying out the integration 
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we obtain after some algebra all the elements in the tensor as 

r ,12 r,i r,i r,i r,i W P V^/ 1 ^ M^Vj + ^0 Oggg^g Mgg + gg^M^ll 

en - [wj - [fcj S 2[fcjB2 - [fcJssFjBs - — 2^ -i J /y, Z k'vo) 2 

r,i r,i T^ ^l^oOgg^M -VoS B3 [k} E1 ) 

612 - - - 2^("1) ^ _ ^ )2 

rH rH r^V S n V k 3( S E3{u] - VoSB2{k]El) 

* = [*w*]«-#»-i)'*^S 



7 - W - fc^o 

r i2 r/i r/i r;i r;i w p W ..^^(-S^M - ^o-SW^si) 
e 22 = [u] - [k] E1 [k] m - [k\ E3 [k\ B3 - — 2_^{-lr -, 77T 

£23 = Wb 2 [A;]b3 

r;l r;l W P W -,\a V oSj 3 S B 2[k] E3 

fl,V 

£32 = [AWfc]B2 

£33 = M 2 - [*M*] B1 - - ^ ^(-^M^L^M 

7 w ^i^o 

(15) 

The dispersion relation is then finally obtained from the condition that 

Det(^) = (16) 

which is valid in any number of dimensions. 

2.2.2. ID and 2D case 

Much can be learned from examining the ID and 2D limits to the general 
dispersion relation. In ID simulations all physical quantities only depend on 
one coordinate Xi, hence [k], k, and k! only have the 1-component. It follows 
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then that the elements of the it are 

,2 



7 



(fa/ - fc^o) 5 



e 22 - e 33 - M - L*J^L*Jsi - — Z.^ 1 ) co'-k\v 



Cl2 — Cl3 — e 21 — ^23 — ^31 — e 32 — 



(17) 



Using Eq. (16), the dispersion relation for the ID case is three uncoupled 
modes, 



en = e 2 2 = 



(18) 



where each mode corresponds to separate components of the electric fields 
Ei, E2, and E3 respectively. Each of these modes is numerically stable as 
long as At is sufficiently small. If we take the limit At — > 0, and As — > 0, 
then Eq. (17) and (18) reduce to the dispersion relations in a real drifting 
plasma (which is completely stable). 



Similarly, the elements of e in the 2D limit can be written as 



oJ 



r, ,12 r,l r.i W V ST ( 1 \H Mggj^ M H + S B3 [g E2^ 

en - M - [k] E2 [k] B 2 - - Ur-1) _ 



ei2 = [^i[A:] B 2-^^(-l) 



M fc^ogji(g g2 M - S B 3Vo[k]El) 
- ^ (u/ - k[v ) 2 



(21 



[k] E 2[k] B1 - - ^(-1) - _ ^- 
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, -r, >1 2 m rjLi -p V7 r ^ S j 2(S E 2[uj] - ggaWgi^o) 
e 22 - M - [fclsiMfl! - - 2^(-l) ^-T7^ 



ri r/i n r;i r;i "P W 1 \m ^(^W^] — ^W^Ui^o) 
e 33 = M - [A;]bi[A;]bi - [^2^52 - 2J -1 ) _ ^ 



ii. i' 



Cl3 — ^23 — ^31 — £32 — 



Using Eq. (16), we can obtain the dispersion relation for the 2D case 

^11^22 — ^12^21 = 633 = 



(19) 



(20) 
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Note that E 3 is de-coupled from the other two directions. 

The numerical features of a particular simulation setup can now be in- 
vestigated by solving the corresponding numerical dispersion relation. Due 
to the use of the finite space and time steps, these dispersion relations not 
only contain terms from the lowest order Brillouin zones (/i = and z7 = 0), 
but also the space aliasing (summation over z7) and time aliasing (summation 
over fi) terms [TJ. The elements of the interpolation dyadic S , and finite 
difference operators [•] comes into the expression due to the finite difference 
treatments when depositing currents and EM fields, and when solving the 
Maxwell Equations. The modifications to the dispersion relation leads to 
numerical instability in the otherwise stable physical system [H 121 [3J H] . 

2. 3. Beam modes and EM modes 

The ID, 2D, and 3D dispersion relations show that a drifting plasma leads 
to beam modes in the dispersion relation that are associated with longitudinal 
oscillations. A beam mode roughly satisfies the dispersion relation 

(u/ - k[v ) 2 = ^~0 (21) 

In addition, a drifting plasma also supports transverse EM waves that 
are described by the dispersion relation 

M 2 - [k] E .{k} B = ^~ (22) 

7 

which for the high gammas are the dispersion relation on the grid for trans- 
verse EM modes in vacuum. In 3D, as the three components of the electric 
field are coupled, we expect to find instabilities near the intersections of the 
EM modes and beam modes in (hi, hi, and ks) space. However, in 2D E3 is 
de-coupled from E\ and E%\ thus the instability (coupling) can only occur in 
k\ and k^ space. 

In ID, all the three components of the electric field are de-coupled so 
instability can only occur if either the longitudinal beam or transverse EM 
mode are numerically unstable themselves. We next concentrate on the 2D 
case since it is easier than the 3D but it still has the possibility of numerically 
unstable modes. 

As we show next, we observe in the simulations and in the solutions to our 
dispersion relation that in fact the unstable wave numbers and frequencies 
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lie at the intersection of u>' — k[vo = (which we refer to as a vacuum beam 
mode) and the vacuum dispersion relation for EM waves. We will use this 
observation to derive asymptotic growth rates for the instability in section |5j 



3. Numerical instability induced by relativistic plasma drift 

3.1. Theoretical analysis of 2D dispersion relation 

Without loss of generality, we use the results in section [2] to study the 
numerical instability induced by the relativistic plasma drift in a 2D system. 
According to the dispersion relation in 2D, we expect to observe instability 
in E 2 (and .83). By calculating the maximum imaginary part of u for real 



values of (ki, k 2 ) for Eq. (20), we can obtain the characteristic pattern of the 
instability in Fourier space, as well as the growth rate of the instability. We 
can also plot the real part of u for ki,k 2 . These results can be used later to 
compare with the simulation results. 

The dispersion relation is general and can be used to examine differ- 
ent choices in Maxwell Equation solvers, in differences between energy and 
momentum conserving field interpolation, in differences between charge con- 
serving and direct current deposition schemes, and the use of smoothing and 
low pass filters. In this paper, we are emphasizing that our dispersion re- 
lation agrees well with the simulation results for cases studied, that we can 
predict the region of unstable modes by plotting where the beam and EM 
modes intersect in k and u space; that we can obtain an asymptotic expres- 
sion for the growth in 3D which agrees well with the simulation for various 
finite difference solvers (including the values of At tht minimize the growth 
rate); and the advantages of using a spectral solver from the point of view of 
eliminating the instability and not attempting to carry out a comprehensive 
survey of all available choices listed above. 

We illustrate the instability using a 2D case with the standard Yee solver 
[15] . We choose the grid parameters and time step that satisfies the Courant 
Condition [2] to eliminate known numerical instability from the EM modes. 
We use the parameters i n Tab. [Tj and substitute the finite difference operator 



for the Yee solver (see Appendix A) into the 2D dispersion relation. We 



assume linear (area) interpolation, momentum conserving field interpolation, 
and a charge conserving current deposition (see Appendix A). 

After obtaining all the roots (u, hi, k 2 ), we plot the dependence of the 
growth rate in the (uj r , ki) space [figure [l] (c)], as well as in the (hi, k 2 ) space 
[figure [I] (d)]. It is evident that all the instabilities are near the main or 
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Parameters Values 

solver Yee 

grid size (fcoAa;i, /C0AX2) (0.1,0.1) 

time step At 0.9xCourant limit 

boundary condition Periodic 

simulation box size (koLi, k^L^) 51.2x25.6 

plasma drifting velocity 7 = 50.0 

plasma density n/n>o = 1 

Table 1: Crucial simulation parameters for the 2D relativistic plasma drift simulation, 

aliased beam modes. Since the terms with < 1, and |^i| < 1 are the most 



important, we neglect higher order terms when solving Eq. (20). Higher 
order /i and v terms can be included in the summation if needed. These ad- 
ditional terms lead to additional unstable modes in (k\, ^2) space with lower 
growth rates as well as to small modifications to the growth rate and location 
of the original modes. A plot in (ki, k 2 ) space with more terms included are 
presented in figure [I] (b), in which we use the asymptotic expression Eq. (30) 
for the growth rate (see section 5] for more details). 

While the results in figure [I (c) and (d) are numerically calculated from 



Eq. (20), the location of the unstable modes can also be conveniently pre- 



dicted by plotting the intersection of the EM modes Eq. (22) and beam 
modes Eq. (21) in (ki, k 2 , co r ) space. This is shown in a 3D plot [figure [T] 
(a)]. By examining the unstable pattern in (ki,k 2 ) space we see that the 
central part of pattern comes from the intersections of the EM modes and 
main beam mode (/i = and v = 0), while the part at the four corners can 
be identified from the intersections of the EM modes and first order spatial 
aliasing beam modes (/z = and ui = ±1). As we argue in section |5j a key 
to mitigating the instability is to manipulate the instability pattern through 
a careful choice of the Maxwell Equation solver. Making a plot in (fc l5 k 2 , 0J r ) 
of the intersection of the EM and beam modes for various solvers becomes a 
useful method for examining where the unstable modes reside without having 
to solve the full dispersion relation. 

3.2. Simulation study of the instability 



To compare with the results in section 3.1, we conducted simulation stud- 
ies in the 2D system using the EM-PIC code OSIRIS [13]. In these simu- 
lations, a neutral plasma with both the ion and electrons drifting in x\ at 
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the same relativistic velocity of 7 = 50.0 is initialized throughout the en- 
tire simulation box. Periodic boundary conditions for fields and particles are 
used. Other crucial parameters for the simulation setup are identical to the 
theoretical study in section |3.1| 

As is shown in figure [2] (a), the total EM energy starts to grow violently 
as the plasma drifts relativistically. The exponential growth indicates that 
a numerical instability is occurring. In addition, the EM field energy in E2 
and B3 and that in E3 and B2 are shown separately. As predicted by the 2D 
dispersion relation the E3 and B2 modes are stable and do not grow. The 
pattern of E 2 at t = 100 u' 1 is plotted in figure fil (e) and (f), and good 
agreement for the location and relative amplitude of the unstable modes is 
obtained when compared against the theoretical prediction [figure [I] (c) and 

(d)]- 

The EM energy grows with a lower rate after t = 110 tu p 1 [figure (a)]. 
The plasma density in this regime is highly modulated by the EM fields. 
The first order perturbation in plasma electron density [figure [2] (b)] shows a 
similar pattern as for E2 [figure [2] (d)] , which confirms they are coupling in 
the system. Note that no exponential energy growth can be seen in the E3 
field [figure g(c)] 

From the simulation we find that for later times after the instability has 
evolved into a nonlinear state, the same pattern in (k%, k 2 ) space as that of the 
linear regime still exists. This indicates that the instability will remain near 
the intersections of the EM modes and beam modes and that both the linear 
and nonlinear growth can be mitigated through eliminating or controlling 
the intersections. 

We also carried out a numerical investigation of the ID dispersion relation 
Eq. (17), and (18) using the same simulation parameters as in Tab. |I] 
(with the ID Courant condition). This confirms that there is no numerical 
instability under these parameters which is expected since E\ is de-coupled 
from E2 and E3 in Eq. (17) and each mode is itself stable. 

We have done the same simulation studies on the use of other finite dif- 
ference solvers beside the Yee solver, including the Karkkainen [HI US] and 
4th-order solver [17]. Good agreement between theory and simulations was 
also found [18] . Some results are shown in figure |5j and discussed in section 
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4. Mitigation of numerical instability 

Due to the fact that the EM dispersion curves for the Yee solver inevitably 
bends down at high k\jk g \ [of. figure [T] (c)] such that the phase velocity of 
electromagnetic (EM) waves on the grid is less than the plasma drifting 
velocity then an instability is found at the low \k\ region in the (A/l,^) 
space due to the fact that the beam mode can intersect the EM modes. In 
addition, there can be instability at high \k\ near the intersection of the EM 
mode and the first order space aliasing beam mode (vi = ±1). According 
to the numerical dispersion relation, the high \k\ part due to aliasing can be 
conveniently smoothed out by applying low pass smoothing to the current 
(and/or EM field) when pushing the particles. However, the part due to 
intersection at the main beam mode is more difficult to mitigate, since the 
physics we want to simulate can reside in this region. 

Results from OSIRIS simulations with the Yee solver, linear order parti- 
cles, and momentum conserving field interpolation with and without current 
smoothing are presented in figure [3] (a)-(c). In the smoothing case 4-pass 
binary smoothing (1,2,1) with 1-pass compensator (-5,14,5) is applied to 
the current. We can see that the part near the first aliasing beam modes 
{yi = ±1) is greatly reduced, while that at the main beam mode {y = 0) in 
the low \k\ region is still present. 

We next explore the use of a spectral solver to eliminate the intersection 
of the EM mode with the main beam mode. The use of spectral solver 
to advance the EM field in Fourier space leads to the EM modes on the 
grid having phase velocities greater than the drifting velocity of the plasma 
which prevents any coupling with the main beam mode [3j [19] . As shown in 
figure [i] (a), the intersection of the EM and beam modes occur first at the 
v\ = ±1 beam mode. Since the EM dispersion surface is a cone in (ki, k2, u) r ) 
plot, we would expect its intersection with the aliasing beam modes to be 
part of an ellipse that resides at high \k\ region. Numerical results obtained 
by solving the corresponding numerical dispersion relation (spectral solver, 
linear shaped particles, momentum conserving field interpolation, and direct 
current deposit) are presented in figure [4] (b) and (c), which confirms the 
empirical prediction in figure [1] (a). 

Moreover, since we are advancing the EM field in the Fourier space, it is 
easy to apply customized filters directly to the EM field. In figure [3] (a), (d), 
(e), we present simulation results using the spectral solver in the UCLA PIC 
Framework [14]. We show results where no filter is used, where a Gaussian 
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shaped low-pass filter is used, and where a low pass filter with a hard cutoff 
is used. The Gaussian shaped filter was of the form, exp(— \k\ 2 /a 2 ), with 
a = 0.9fci the RMS width of the filter in Fourier space, and the low pass 
filter with a hard cutoff set modes with \k\ > 0.9fci to zero. As we can see in 
the hard cutoff case, the instability at V\ = ±1 is essentially eliminated from 
the simulation, thus only leaving instability with higher order terms. This is 
confirmed by the fact that the energy of the transverse EM modes remains 
at a very low level [figure [3] (a)]. 

5. Asymptotic expression for instability growth rate 

5.1. Derivation of asymptotic expression 

In section |3j we obtained the instability pattern and growth rate by nu- 
merically solving the dispersion relation equation, which is feasible on a mod- 
ern laptop computer in 2D scenario, but much more difficult in 3D. As men- 
tioned above, we observed the highest growth rate at the intersections of 
the beam modes and EM modes. Taking advantage of this observation, we 
here derive an asymptotic expression for the solutions near the beam mode 
oj' = k[v . These expressions will not only speed up the instability pat- 
tern analysis in 3D, but also provide more insights into the dependence of 
instability pattern and growth rate to the grid sizes and time step used in 
simulation. 

We denote the terms which are summing over fi and V in as Qij, i.e. 

en = [u] 2 - [k] E2 [k] B2 - \k\m\k\m - Qu e i2 = [k] E i[k] B 2 - Q12 e 13 = [k] E1 [k} B3 - Q 13 
£21 = [k] E 2[k)m - Q21 £22 = [uf - [k]Ei[k]Bi ~ [k} E 3[k]m - Q22 £23 = [k] E 2[k]m 
£31 = [k] E 3[k] B1 - Q31 £32 = [k} E 3[k] B2 £33 = M 2 - [k] E1 [k] B1 - [k] E2 [k] B2 - Q33 

We expand u' around the beam mode u>' = k[v , and write co' = k[v + 5u', 
where 5u' is a small term. In addition, we will use the relativistic limit 
Vq — y 1. We will also truncate the summation over v 2 and ^3, keeping only 
the v 2 = V3 = terms. Using det( e ) = 0, and dropping terms of higher 
order of (u 2 /^) 2 , (u 2 /-y) 3 , . . ., we can obtain 

A x bJ 2 + B x 5oj' + d = (23) 
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where 



A 1 = [u] 2 ([uf - [k] E1 [k]m - [k] E2 [k] B2 - [k) E3 [k) B3 ) 

Bi = [k] E1 [k} B2 Q 21 + [k] El [k] B3 Q 31 - ([u] 2 - [k\ E2 [k\ B2 ) Q 22 - ([tu] 2 - [k] E3 [k] B3 ) Q 33 
d = - (H 2 - [k] E1 [k] m ) Q u + [A^ 2 [A;] B1 Q 12 + [k) E3 [k) B1 Q 13 (24) 

Now we use the condition that (u/, k[) sits near the EM modes, 

N 2 « + [A;] B2 [A;] B2 + [^[A^a (25) 

We further expand the finite difference operator [cj] = £ + Soj'^i, where 

sin(A;At/2) ,~ A . . 

e ° = At/2 6 = C0S (^At/2) (26) 

and 

k = k\ + — /iCJg (27) 

We further expand [w] to first order in Ai since this term is sensitive near 
the EM mode, while keeping only zero order of [w] in B 1: and C\. We then 
obtain a cubic equation 

A 2 5uj' 3 + B 2 5u' 2 + C 2 Su' + D 2 = (28) 

with the coefficients 
A2 = 2£ 3 6 

B2 = e fe 2 - ([k}Ei[k]m + [k] E2 [k} B2 + [k} E3 [k] B3 )} 

C 2 = [k] E1 [k] B2 Q 21 + [k] E1 [k) B3 Q 31 - - [k) E2 [k) B2 ) Q 22 - (e - [k] E3 [k] B3 ) Q33 
D 2 = - - [k] E1 [k]m) Q u + [k} E2 [k} m Q 12 + [k] E3 [k] B1 Q 13 (29) 

When calculating the instability growth rate, we obtain the imaginary part 
of roots ^s{5u'(k, fi, for each /1 and v\ by solving Eq. (j28~])-(29 ), and the 
growth rate r(/c ) for a particular mode k is chosen to be max{^s{Su'(k , /i, z^i)}}. 
When solving for each ^s{Suj'(k, we only keep the corresponding /i and 

z/i terms in the above cubic equation. Eq. (28)-(29) can be used to plot the 
growth rate in Fourier space, and can be conveniently simplified to 2D. In fig- 
ure [T](b) we plot the asymptotic instability growth rate with /i = 0, \v±\ < 4 
for 2D Yee solver using the same parameter listed in table [TJ 
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Near the transverse EM modes [u>] 2 ~ £q ~ [k] E ■ [k] B , we can drop the B 2 



term in Eq. (28). According to our numerical results, we can further simplify 



the analytical expressions by dropping the small C 2 term. The asymptotic 
growth rate T(k) for this mode corresponds to the maximum imaginary part 
for the three roots, 



r(*) 



(Co - [k]Ei[k]m)Qu ~ [k]E2[k] B iQvi - [^W^siQia 



1/3 



V3 • 2 
2 



u*Sji {(S B3 £o ~ S E 2[k] B1 )[k} E2 k2 + (S B2 £ - S E3 [k] B1 )[k} E3 k 3 } 



27fo 2 fi 



1/3 



(30) 



This expression shows the relation between the instability growth rate and 
grid sizes, time step, interpolation and smoothing functions, and finite dif- 
ference operators. Note that from the positions of the interpolation func- 
tions we could immediately see that using higher particle shape, or using a 
stronger smoother helps mitigate the instability pattern at high \k\ region, 
which agrees well with our simulation. 

5.2. Parameter scans for minimal instability growth rate 



With the asymptotic expression Eq. (28)-(29), and also Eq. (30), we can 



greatly speed up the solution of numerical dispersion relation in 2D and 3D. 
In addition, the asymptotic expression makes the parameter scan to study 
the dependence of instability pattern and growth rate between various grid 
sizes and time step more convenient. 

In figure|5j we scanned the grid sizes Ax\ and time step At/ Ax\ for the 2D 
and 3D Yee solver, and Karkkainen solver, and compared the growth rates 
with the OSIRIS simulations. We have kept Ax± = Ax 2 {= Ax 3 ) during 
the parameter scan for 2D (and 3D). We likewise plotted out the OSIRIS 
simulation data for Ax\ = 0.1 together with the asymptotic data (similar to 
the plots in Ref. [S])- There are several interesting points worth noting in 
figure [5] First, we can see there is a "magic time step" [5] At m /Axi where 
the growth rate is minimized in most cases; on the other hand, the instability 
growth rate decreases monotonically as the grid sizes increases; second, when 
the grid sizes are square (2D) or cubic (3D), the "magic time step" At m /Axi 
is an invariant for different Axi, in both the momentum conserving (MC) 
scheme, and energy conserving (EC) scheme; third, the instability growth 
rate for 2D and 3D are nearly the same for given Ax± and At/Ax\ under 
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the same field interpolation scheme; the values for the magic time steps are 
also nearly the same in 2D and 3D (note that according to the asymptotic 
expression, the magic time step for Yee solver 3D EC scheme also resides at 
around At m /Axi ~ 0.65, but we did not plot it out since that At m is beyond 
the Courant limit for this solver). The parameter scan using the asymptotic 
expression for the Karkkainen solver with EC scheme shows the "magic time 
step" at around At m /Axi = 0.7, which agrees with the results reported in 
Ref. [8]. However, according to our simulation and theoretical results, we 
found the "magic time step" not only in Karkkainen solver, but also in Yee 
solver; and not only for EC scheme, but also for MC scheme. This is also 
reported in [6] for the 2D cases. 

The fact that the "magic time step" At m /Axi does not depend on Ax\ 
for square (and cubic) cell for the Yee and Karkkainen solver is evident from 
Eq. ( |30~) ). Applying the detailed form of finite difference operators [k]Bi for 
Yee and Karkkainen solver (note they have the same [fcjsi, see Appendix A[ ), 



for both MC and EC scheme, the expression of T can be expressed as lA^) 1 ' 3 
times a function of At/Ax\, and k^/kgi. (this function is different for different 
field interpolation schemes). Since k'Jkgi ranges from (—0.5,0.5) regardless 
of Axi when calculating the growth rate, the extreme value of T resides at 
the same At/Axi for different Axi (although different field interpolation 
schemes give different "magic time step"). In particular, in MC scheme the 
terms 

Sb3^,o ~ SW^]m Sb2^,o — SW^]m (31) 



or only the first one in 2D) in Eq. (30) are zero when At/Ax\ = 1/2 for 



these two solvers. As a result, both Yee and Karkkainen solver reach the 
minimal growth rate at At/Ax\ = 1/2 in MC scheme, which agrees well 
with OSIRIS simulations. 



6. Conclusions and future work 

We derived a general multi-dimensional numerical dispersion relation for 
the relativistic plasma drift in the EM-PIC simulation that can include dif- 
ferent choices in Maxwell solvers, differences between energy and momentum 
conserving field interpolation, differences between charge conserving and di- 
rect current deposition schemes, and the use of smoothing and low pass 
filters. In this paper we emphasized trying to understand the source of the 
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instability and the structures of the dispersion relation. We confirmed that 
no instability occurs in ID, and that in 2D and 3D a strong instability occurs 
due to the coupling between beam modes and transverse EM modes in the 
system. We can predict the pattern and growth rate of the instability for a 
particular simulation by solving the corresponding numerical dispersion rela- 
tions. An asymptotic expression which permits rapid parameter scan for the 
ranges of unstable modes was derives. These results are compared against 
simulation results using the EM-PIC code OSIRIS [13], as well as UCLA PIC 
framework [14], and good agreement is obtained. 

Moreover, by plotting the intersection of the EM and beam modes in 
(ki, k 2 ,w r ) space the numerical instability patterns can be conveniently pre- 
dicted. Maxwell Solvers, such as spectral solvers, that have EM waves with 
phase velocities greater than the plasma drifting velocity are therefore free 
of any instability due to the main beam mode. We showed that the use of 
a spectral solver does indeed eliminate any instability from the main beam 
mode, leading to instability predominantly from coupling between the EM 
wave and the lowest order aliased beam modes. This coupling exists only 
at high | A; | area at the edge of the lowest order Brillouin zone. We showed 
that a low pass filter with a hard cutoff can eliminate these modes without 
effecting lower \k\ modes that are physically important in a properly resolved 
simulation. 

In addition, by using the fact that the modes with highest growth rate are 
found near the intersection of the beam modes and EM modes, we derived 
an asymptotic expression for the growth rate that can be useful to study 
the growth rates with various smoothing functions and different cell sizes. 
The asymptotic expression speeds up the calculation of instability growth 
rate and makes the investigation of instability pattern and growth rate in 3D 
feasible. By conducting parameters scan using the asymptotic expression for 
the Yee, and Karkkainen solver in 2D and 3D, we confirmed the "magic time 
step" that minimize the growth rate, as reported in [8j. We found the ratio of 
the magic time step over the grid size along the drifting direction, At m /Axi 
is determined by the field interpolation scheme used in the simulation, yet 
stays constant for different Ax\. These observations agree well with the 
simulation results. 

This paper reports on efforts to understand and mitigate the instability 
in cases in which there is only a drifting plasma. Areas for future work 
include trying to optimize the particle order, field smoother, and field solver 
for mitigating this numerical instability in LWFA simulation, understanding 
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how the various solvers and use of specific time steps effect the numerical 
dispersion relation for light waves, and for studying if the use of the Yee 
mesh together with a spectral solver also leads to an optimal time step. 
In addition, future work should include exploring the tradeoffs in accuracy, 
computational speed, and parallel scalability for the different choices. 

This work was supported by DOE awards DE-FC02-07ER41500, |de-sc0008491 
DE-FG02-92ER40727, and [de-sc0008316[ and by NSF grants NSF PHY- 
0904039 and NSF PHY-0936266. Simulations were carried out on the UCLA 
Hoffman 2 Cluster. 

Appendix A. Interpolation dyadic and finite difference operator 

In this appendix we will write out the explicit expressions for the in- 
terpolation dyadics o for the fields and the currents. For a momentum 
conserving scheme, in 3D the interpolation dyadic for the EM field after the 
Fourier transform can be expressed as: 

Sex = si t isi t 2Si }3 r)i S E 2 = si tl si j2 si t3 T]2 S E3 = si tl si j2 si >3 r] 3 

S m = COs(wAt/2)s; 5 iS />2 S Zi 3?72?73 S B 2 = COs(u At /2) S Ij iSj, 2 Sj ) 3»7i773 

S B3 = cos(wAt /2)si t isi^si t3 r] 1 r] 2 (A.l) 



and rji = ( Ui } ( = — 1 when the EM field has half-grid offset in the i direction, 
and ( = 1 when it is defined at grid point (this was the term missing in our 
earlier version). / refers to the order (I = 1 is area weighting or linear 
interpolation for the charge). 

For an energy conserving scheme, we have 

S E 1 = Si-i t iSi t 2Si j3 T]i S E 2 = Si t iSi- lj2 Si t 3r)2 S E 3 = Sj^Sj^Si-l^^ 

S B i = cos(u At/ '2) si A si_ lt 2Si^ 3 r] 2 r]3 S B2 = cos(uAt/2)si_ 1:1 si t2 si- 1:3 r] 1 r] 3 
S B3 = cos(u}At/2)si- 1A si- li2 si i3 ri 1 r)2 (A. 3) 

For the rigorous charge conserving scheme (as is used in OSIRIS), the current 
interpolation dyadic is approximately: 

S n = Si-i tl Si )2 Si !3 T)i S j2 = S Zi iS;_i )2 Sj i3 ?72 ^3 = S^iSj^Si-l^ ( A - 4 ) 



where 




(A.2) 
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We note that the expressions for Sj are strictly valid for the charge conserving 
scheme in the limit of At — > 0. Meanwhile, when the current is directly 
deposited (as is done in the UPIC framework), the current interpolation 
functions are, 

Sji = si^s^Si^rjx S j2 = si A si t2 si t3 r] 2 S j3 = Sj^s^sj^ (A.5) 
The space finite difference operator for the Yee solver is: 

[k]j = ^2) (A 6) 

and is the same for electric and magnetic field. The space finite difference 
operator for the Karkkainen solver is the same as Yee solver for the magnetic 
field, as for the electric field, we used 

sin(fc i Ax i /2) 
[k] Ei = Ci A ^ /2 (A.7) 

where 

c i = Qi + 2# 2 {cos(fc 2 Ax 2 ) + cos(fc 3 Aa; 3 )} + 4# 3 cos(fc 2 A:r 2 ) cos(/c 3 Ax 3 ) 

c 2 = 0i + 2^ 2 {cos(/c 3 Ax 3 ) + cos(fciAa;i)} + 4# 3 cos(fc 3 Ao: 3 ) cos(/ciAa;i) 

c 3 = Q\ + 26 l 2 {cos(A;iAxi) + cos(fc 2 Aa; 2 )} + 4# 3 cos(fciAxi) cos(/c 2 Aa; 2 ) 

(A.8) 

and 

6 l = 7/12 9 2 = 1/12 fl 3 = 1/48 (A.9) 

are the tunable parameters for the Karkkainen solver |16| . The space finite 
difference operator for the spectral solver is 

[k]i = h (A. 10) 

The time finite difference operator for the Yee, Karkkainen, and spectral 
solver are the same 
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Figure 1: Numerical instability pattern in the Yee solver. Growth rate are color-coded, 
and normalized with w g . (a) EM modes intersect with the main beam mode (/i = 0, 
v = 0), and first order space aliasing modes (/i = 0, v\ = ±1); (b) is the instability 
pattern = 0, |^i| < 4) in (fci,fc 2 ) space, plotted using Eq. p8]H29l; (c) and (d) are 
the instability pattern < 1, \ui\ < 1) in (ui r ,ki) and (fci,fc 2 ) spaces obtained from 



solving Eq. (19) and (20). EM modes for different propagating angles [in degree] and the 
beam modes are likewise plotted in (c). (e)2^esents the corresponding simulation results 
in (uj r , fci) space, and (f) in (fei, ^2) space. Data in (e) and (f) show the modes present at 
t = 100 u)~ , and are not a measurement of the growth rates. 



(a) particle(blue); E3/B2 (green); E2/B3 (red) ( ' electron density perturbation 




k /k j k 7k . 

1 gi 1 gi 

Figure 2: We present in (a) the energy evolution of the EM energy for the two cases. 
The corresponding dotted line indicates their variation in time after t — 100 w" 1 ; (b) is 
the plasma electron density perturbation in (fci,A;2) space, (c) presents the E% in {k\,k%) 
space, and (d) presents the Ei in (ki,k-z) space. 
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(a) 




k./k , k„/k . 

1 g1 1 g1 

Figure 3: Instability mitigation in 2D simulation for the Yee and spectral solver, (a) 
Energy evolution of the E2 field in various simulation setups; (b) presents E2 in (fci,/^) 
space for the Yee solver case in which a 4-pass binary smoothing with compensator is 
applied to the current, while (d) is the E2 in (ki, fe) space for spectral solver with cutoff 
smoothing, (c) and (e) is the E2 field in (kx, £2) space for the non-smoothing case for the 
Yee, and spectral solver respectively, (b)-(e) are plotted at t = 240 u~ x . 
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Figure 4: Numerical instability pattern for the spectral solver. Growth rate are color- 
coded, and normalized with u) g . (a) EM mode intersects with the first order space aliasing 
mode (p = 0, V\ = ±1); (b) and (c) is the instability pattern in (bj r ,ki), and (fci,^) 
spaces obtained from solving Eq. (19) and (20). Beam modes (fj, = 0, |^i| < 1) and EM 
modes for different propagating angles [in degree] are likewise plotted in (b). 
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Figure 5: Parameter scan of Ax\ and At/ Ax\ for the Yee (first two rows) , and Karkkainen 
(last two rows) solver. The first and third row uses momentum conserving (MC) scheme, 
while the second and fourth row uses energy conserving (EC) scheme. The simulation 
results are likewise plotted in (c), (f), (i), and (1) at Axi — 0.1 for comparisons. In (c) 
and (f) the dotted line at At/ Ax\ w 0.577 is the 3D Courant limit (CL), and that at at 
At/Ax-i « 0.707 is the 2D CL. 
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